地図データ例

library(maptools)
map <- readShapePoly(system.file("etc/shapes/columbus.shp",package="spdep")[1])
plot(map)

複数のポリゴンの融合

rgeosパッケージのgUnionCascaded関数を用いて複数のポリゴンを融合して、1つのポリゴンを作成する。

library(rgeos)
map.u <- gUnionCascaded(map) # 融合したポリゴンを出力
plot(map.u)

2つのポリゴンの共通部分

rgeosパッケージのgIntersection関数を用いて、2つのポリゴンの共通部分のポリゴンを作成する。 例として、下図の赤線で描いた四角領域と、地図ポリゴンの共通部分を赤色で塗りつぶす。

# 例の四角領域ポリゴンの作成
bb <- bbox(map.u)
bb[,2] <- bb[,2]*0.8
box <- Polygon(cbind(c(bb[1,1:2],bb[1,2:1],bb[1,1]),
                     c(bb[2,1],bb[2,1:2],bb[2,2:1])))
box <- Polygons(list(box),ID="box")
box <- SpatialPolygons(list(box))
plot(map.u)
plot(box,add=TRUE,border="red")

plot(map.u)
int <- gIntersection(map.u,box) # 共通部分のポリゴンを出力
plot(int,col="red",add=TRUE)